# 符号计算带电球的电势能
# 注释部分由AI完成
# \documentclass[12pt, a4paper]{article}
# \usepackage{geometry}
# \geometry{
# 	a4paper,
# 	left=12.7 mm,
# 	right=12.7 mm,
# 	top=12.7 mm,
# 	bottom=12.7 mm,
# }
# \usepackage{graphicx}
# \usepackage{amsmath}
# \usepackage{physics}
# \graphicspath{ {./} }
# \usepackage{ctex}
# \usepackage{verbatim}
# \usepackage{makecell}
# 
# \newcommand{\bvec}[1]{\ensuremath{\mathbf{#1}}}
# 
# \begin{document}
# 考虑一个半径为 $R$、总电荷为 $Q$ 的均匀带电球体，计算其静电能。设介电常数为 $\epsilon$（代码中为 \texttt{epi}）。
# \subsection*{电场分布}
# 根据电场的Gauss定理
# \begin{equation}
# 	\div \bvec E = \frac{\rho}{\epsilon}
# \end{equation}
# 两边同时对半径为$r$的球体积分
# \begin{equation}
# 	\int \div \bvec E \dd V = \int \frac{\rho}{\epsilon} \dd V = \frac{Q_{in}}{\epsilon}
# \end{equation}
# 根据散度定理，
# \begin{equation}
#  \int \bvec E \cdot \dd \bvec A  = \frac{Q_{in}}{\epsilon}
# \end{equation}
# 考虑到$E$的球对称性
# \begin{equation}
# 	E \cdot (4 \pi r^2)  = \frac{Q_{in}}{\epsilon}
# 	\Rightarrow E = \frac{Q_{in}}{4 \pi r^2 \epsilon}
# \end{equation}
# 随后需要分类讨论$Q_{in}$。
# \begin{equation}
# 	Q_{in} = 
# 	\begin{cases}
# 		Q, r > R\\
# 		\frac{\frac{4}{3} \pi r^3}{\frac{4}{3} \pi R^3} Q = \frac{r^3}{R^3} Q, r<R\\
# 	\end{cases}
# \end{equation}
# 因此，
# 球内 ($r \leq R$) 和球外 ($r > R$) 的电场分别为：
# \begin{equation}
# 	E_{\text{in}} = \frac{Qr}{4\pi\epsilon R^3}, \quad E_{\text{out}} = \frac{Q}{4\pi\epsilon r^2}
# \end{equation}
# \subsection*{能量密度}
# 静电能密度由下式给出：
# \begin{equation}
# 	u = \frac{\epsilon}{2}E^2
# \end{equation}
# 因此球内外的能量密度为：
# \begin{equation}
# 	u_{\text{in}} = \frac{\epsilon}{2}\left(\frac{Qr}{4\pi\epsilon R^3}\right)^2, \quad u_{\text{out}} = \frac{\epsilon}{2}\left(\frac{Q}{4\pi\epsilon r^2}\right)^2
# \end{equation}
# \subsection*{能量计算}
# 总静电能通过对能量密度的体积分得到，此处使用球坐标积分：
# \begin{equation}
# 	U = \iiint (u_{in} + u_{out}) dV = U_{\text{in}} + U_{\text{out}}
# \end{equation}
# 球内部分 ($r \leq R$)：
# \begin{equation}
# 	U_{\text{in}} = \int_0^{2\pi} \int_0^\pi \int_0^R u_{\text{in}} r^2 \sin\theta dr d\theta d\phi = \frac{Q^2}{40\pi\epsilon R}
# \end{equation}
# 球外部分 ($r > R$)：
# \begin{equation}
# 	U_{\text{out}} = \int_0^{2\pi} \int_0^\pi \int_R^\infty u_{\text{out}} r^2 \sin\theta dr d\theta d\phi = \frac{Q^2}{8\pi\epsilon R}
# \end{equation}
# 总静电能：
# \begin{equation}
# 	U = U_{\text{in}} + U_{\text{out}} = \frac{3Q^2}{20\pi\epsilon R}
# \end{equation}
# \end{document}
# Gitee Repo

import sympy

Q,r,R,epi = sympy.symbols('Q,r,R,epi')
theta,phi = sympy.symbols('theta,phi')

# 内部
Ein = Q*r/(4*sympy.pi*epi*R**3)
uin = epi/2*Ein**2

# 外部
Eout = Q/(4*sympy.pi*epi*r**2)
uout = epi/2*Eout**2

# 球积分
Uin = sympy.integrate(uin*r**2*sympy.sin(theta),(r,0,R),(theta,0,sympy.pi),(phi,0,2*sympy.pi))
Uout = sympy.integrate(uout*r**2*sympy.sin(theta),(r,R,sympy.oo),(theta,0,sympy.pi),(phi,0,2*sympy.pi))
U = Uin+Uout

print(Uin,Uout,U)

# 结果 :
# U_in  = Q**2/(40*pi*R*epi),
# U_out = Q**2/(8*pi*R*epi),
# U_tot = 3*Q**2/(20*pi*R*epi)
